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The pin-loaded lug is a structural element of considerable importance in 
aircraft design, particularly in the design of helicopter rotor and control 
systems. Much work has been done in the analysis of lugs subjected to 
static loads. As a result, the static analysis of lugs has been reduced to 
a well-established rational convention, the most notable work being a much- 
referenced paper by Melcon and Hoblit wherein design allowables and an 
interaction formila for statically loaded aluminum and steel alloy lugs are 
reported. In contrast to the static case, no analogous design criterion 
exists for the design of lugs simultaneously subjected to axial and trans- 
verse fatigue loads. A most glaring testimony to the deartl of valid experi- 
mental data on pin-loaded lugs is demonstrated in MIL-HDBK-5A, wherein the 
section on joints offers no design guidance for lugs. 


This contract was initiated to: 


@ Evaluate the fatigue strength of lugs subjected to vibratory loacings 
at various orientations to the lug axis of symmetry. More specifically, 
an interaction formula relating load orientation to lug endurance Limit 
was sought. 


@ Substantiate the photoelastically established benefits in lug fatigue 
strength that can be derived through selection of interference fit. 


@ Determine the influence of edge distance and material on lug fatigue 
strength. 


Seventy-three lug specimens were validly failed by step-testing leading to 
the development of design charts in the corm of modified Goodman diagrams 
for each material, load direction, and interference fit at two probability- 
of-failure levels. These charts compare favorably with test results 
reported in the literature and satisfy structural requirements for a range 
of edge-distance and load ratios particularly suitable for use in helicopter 
design. Development of an interaction formula did not materialize. Ex- 
cessive scatter in the data precluded development of a general interaction 
formula applicable to both steel and titanium for each edge-distance ratio 
tested. A specific interaction formula for each contiguration tested, 
although possible, was not pursued. 


Results conclusively demonstrate that lug fatigue strength is materially 
improved by the introduction of high interference fit. Verification of the 
existence of an optimum interference fit as photoelastically predicted was 
inconclusive. For the high-modulus mterials tested, the level of inter- 
ference obtainable was limited by attainable thermal size changes. Thus, 
the “optimum was the maximm attainable interference fit nvt causing lug 
yield. For lugs of lower modulus, such as aluminum or steel and titanium 
lugs with liners having substantially heavier wall thickness, it is believed 
that an optimum interference fit does exist beyond which increased inter- 
ference would be detrimental. 


Task 1F162204A14601 
Contract DAAJ02-67-C-0066 
USAAVLABS Technical Report 70-49B 
November 1970 


FATIGUE STRENGTH OF LUGS 
CONTAINING LINERS 


VOLUME II 
COMPUTER PROGRAM 
USED FOR ANALYSES 


By 
Robert Mayerjak 


Prepared by 


Kaman Aerospace Corporation 
Bloomfield, Connecticut 


For 


U. S. ARMY AVIATION MATERIEL LABORATORIES 
FORT EUSTIS, VIRGINIA 


This document is subject to special export controls, and each 
transmittal to foreign governments or foreign nationals, may 
be made only with prior approval of U. S. Army Aviation 

Materiel Laboratories, Fort Eustis, Virginia 23604. 


SUMMARY 


This report presents a FORTRAN program for the analysis of 
elastic, two-dimensional, plane-stress structures, Examples 
show the application of the program to the analysis of lugs. 


iii 


FOREWORD 


This project was performed under Contract DAAJ02-67-C~-0066, 
Task 1F162204A14601, under the cognizance of Mr. Joseph H. 
McGarvey of the Aeromechanics Division of USAAVLABS. 


The tests and analyses were conducted at the Kaman Aerospace 
Corporation, 


The report consists of two volumes: 
Volume I, Results 


Volume II, Computer Program Used for 
Analyses 


The computer program presented herein was developed with 
contractor funds prior to the contract. The very significant 
contributions of Mr. Alex Berman and Dr. John Hsu to the 
computer program are gratefully acknowledged, 
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INTRODUCTION 


This volume presents a FORTRAN program which was used to 
calculate the stresses in the lugs reported in Volume I. 

The program is named MA2B, which is an abbreviation for matrix 
analysis, two-dimensional structures, program version B. The 
information reported herein will enable the reader to use the 
program and to modify it, if desired. 


MA2B is itself a modification of another program*. This 
reference provides excellent, lucid derivations and descrip- 
tions of the structural analysis methods and the basic logic 
used in programming. It is believed that it contained the 
first well-documented, general-purpose structural analysis 
program made available without proprietary restrictions. This 
volume draws heavily from the referenced document. 


Program MA2B differs from its parent in four ways: 
1. It contains a general, quadrilateral-shaped element. 


2. It uses a special Gauss-elimination algorithm for 
the solution of the simultaneous equations. 


3. It is restricted to two-dimensional analysis. 
4. It uses more compact input and output. 


The incorporation of a general, quadrilateral-shaped element 
is particularly important for problems with curved boundaries, 
such as lugs. For these problems the quadrilateral element 
gives more accurate and more easily interpreted results. The 
special Gauss-elimination algorithm greatly increased the 
number of nodes which could be used for a given core memory 
size. The algorithm stores only the non-zero elements of 

the upper triangle of the master stiffness matrix plus up to 
five columns of non-zero forces. Two facts are noteworthy 
regarding this technique: 


* Przemieniecki, J.S., and Berke, Laszlo, DIGITAL COMPUTER 
PROGRAM FOR THE ANALYSIS OF AEROSPACE STRUCTURES BY THE 
MATRIX DISPLACEMENT METHOD, FDL TDR 64-18, AF Flight 
Dynamics Laboratory, Wright-Patterson Air Force Base, Ohio, 
April 1964, AD600418. 


1. The economy of storage exists even if a few elements 
in the stiffness matrix are very widely spaced (not 
banded about the diagonal). 


2. The number of nodes that can be used increases in 
linear proportion to the size of the available core. 


MA2B was restricted to two-dimensional analysis to increase 
further the number of nodes that could be used, As a result 
of the changes, MA2B can solve problems with 300 node points 
using double precision arithmetic. For comparison, the 
referenced program allowed only 70 nodes and used single 
precision arithmetic. As a consequence of the increased size 
of the problem that could be handled, it became desirable to 
Goose the input and output to more compact block tabular 
orms. 


STIFFNESS MATRICES FOR QUADRILATERAL ELEMENT 


The referenced document provides a complete presentation of 


the method used, including derivations of all equations. This 


section describes only the additions that were made. 


The coordinate system used for the quadrilateral element is 
shown in Figure l. 


Figure 1. Coordinates for Quadrilateral Element. 


The stiffness matrix for the quadrilateral element was 
calculated using the essumed displacement function: 
uy= D, + Dad + Day + Dy ty 
Any © Dd, + Dea + Diy + Dety 


If P-Q is not perpendicular to R-S, the constants can be 
determined from the given coordinates for points P-Q-R-S, 
called 1-2-3-4, respectively. Thus, 
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It can be shown that all elements of F are zero except 


vs (+545 — Wels. * 4a4s(te a Oy ye 
fa les * 4344 (ts- 4a)/ Cp 42) 

Fact ta yep 
nals = — bays /p 
F,=-Fag* G2 eh = - O/rga 
face ge? (Cx 43) a Cape ey 
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where = 437% (rga-ya) and the hierarchy of 
operations of FORTRAN apply. 


The strains become 
Cry O Aby D, + D, ay 
€ = Cast = Ouy = D, + Dat 


Cny omg 1 ity Ds + Dee +D, +Dey 


e= Au 


(, + Fs) (E+ Fay) aoe (Rete) 
A= (e. + Ar) (F.+ 82 4) aes (Ge + Fee 4) 
(Fs; +f, z Cd + Ferg)" (Esthg+ ya + bey 


The element stiffness matrices k and h were calculated using 
the unit displacement theoren. 


v 
The matrix A was considered to be composed of 3 components. 


A =A, +A,n+Asy 
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INSTRUCTIONS FOR USER 


This section provides instructions for using the program and 
descvibes in detail the input and output. Complete examples 
are given to demonstrate the use of the progran. 


GENERAL 


The following conments provide general guidance for input 
preparation. 


Units 
Any consistent set of units may be used for the input. The 


output will have corresponding units. The following units 
are recommended: 


input: Loads, kips output: forces, kips 
lengths, in. deflections, in. 
E, ksi stresses, ksi 


Shape of Elements 


The program will run to completion with almost any combination 
of elements. The only element which kills the program is a 
special quadrilateral element described in more detail later 
in these comments. In most cases, the results will have 
engineering significance even for very oddly shaped elements. 
However, the best results for stress analysis will be 
obtained if the following rules are applied: 


1. Avoid panel elements with large length-to-width 
ratios. : 

2. Avoid mixing triangular and quadrilateral elements 
in the region of interest. The triangular elements 
will appear to be stiffer and will disturb the stress 


distribution. 

3. Use a gradual transition from large-to small-sized 
elements. 

4. Use quadrilateral elements which are nearly 
rectangles. 


5. Orient the local axes of the quadrilateral elements 
in the direction of anticipated principal stresses. 


Nusber of Nodes 


The permissible number of nodes is limited by the number of 
non-zero elements that are developed in the stiffnese matrix 
and the loading conditions, both initially and Gareat the 
course of the Gaussian elimination. The program is dimen- 
sioned to handle 6000 non-zero elements. Now many nodee thie 
corresponds to is indeterminate. However, for a typical 
structure, 300 node points can be used if the nodes are 
numbered judiciously. For example, the lug described herein 
had 258 node points and developed less than 4200 non-sero 
elements. The objective in node-point numbering is to avoid 
interspersing high and low node-point numbers. A 

practice is to assign the node numbers consecutively, starting 
at one end of the structure and proceeding systematically 

to the other end. This will cause the stiffness matrix to 
become a desirable narrow band along the diagona. and will 
make the results more easily interpreted. If the structure 
is closed, the above-described procedure will cause the 
highest node numbers to become adjacent to the lowest numbers. 
This is not to be feared. It will not appreciably disturb 
the Gaussian elimination procedure used herein. 


Number of Elements 
Any number of elements can be used. 


Number of Loading Cases 


Up to five peony Senne can be analyzed simultaneously, with 
a small increase computation time. However, it is impor- 
tant to note that a large number of non-zero terms can be 
created if several loading conditions are used which have 
many node points loaded mechanically or by thermal forces. 
Por such cases, the permissible AG@aber of nodes ie lees than 
if the load cases were done one at a time. 


Datum Coordinates 


Any convenient origin for datum coordinates may be used. 
Either right-hand or left-hand systems may be used. 


Support Conditions 


The support conditions are specified by placing the letter 5& 
after the coordinate, in a proper columm. The computer wil) 
then prohibit the displacement of that node point in the 


direction of the coordinate. Sufficient supports must be 
defined to enable the structure to be stable for any loading 
condition. 


Mechanical Loads 


The actual mechanical loads acting upon the real structure 
must be applied as concentrated forces acting at the node 
points of the idealized structure. The directions of the 
forces must correspond to datum coordinates, not local 
coordinates. 


Thecmal Loads 


Thermal stresses can be determined by specifying a coefficient 
of thermal expansion and a temperature for each element. The 
tesperature distribution in the real structure must be 
represented by temperatures which can vary from element to 
Clement, but which are constant within each element. If no 
thermal loads are applied, simply leave the input spaces 
blank. 


Element Identification 


The element identification number IE may be assigned randomly. 
However, a systematic assignment will be very helpful for 
finding the element during the interpretation of the output. 


Local Coordinates 


The node numbers of the vertices of each element are used to 
Gefine both the position of the element and a set of local 
coordinates. The node numbers appear in the input as IP, IQ, 
IR, and I8; they correspond to the points P-Q-R-S shown in 
Pigure 2. 


Figure 2. Local Coordinates. 


The local coordinates xy shown in Figure 2 conform to the 
following rules: For bar elements, local x has the direction 
and sense of a line drawn from P to Q. For panel elements, 

y has the direction and sense of a line drawn from P to Q. 
Local x is perpendicular to y. Local x has a sense that can 
‘be established using the right-hand rule, considering local 

z to be coming out of the sheet. 


The direction of local coordinates is important because the 
stresses will be calculated and presented in the output using 
the local coordinate system. 


A_Special, Requirement for a Quadréiateral 


P-Q-R-S for a quadrilateral must be selected in consecutive 
order around the element. Best results will be obtained if 
R-S is relatively parallel to P-Q. If the R-S direction were 
chosen to be perpendicular to the P-Q direction, as shown in 
Figure 3, the method of analysis would break down. Under 
such conditions, the coefficients D for the assumed displace- 
ment function could not be determined uniquely from the 
coordinates of the vertices of the quadrilateral element. 

The matrix F would be singular. The condition can always be 
avoided by a simple change of designation for P-Q-R-S, as 
shown in Figure 3. 


Very Bad, Satisfactory 
Fatal to Program 


Figure 3. Node Identifications P-Q-R-S 
for Quadrilateral Element. 
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OUTPUT DESCRIPTION 


Name Description 

XX normal stress parallel to the local x~-coordinate 

YY normal stress parallel to the local y-coordinate 

XY shear stress 

ON octahedral normal stress 

Os octahedral shear stress 

NZE non-zero elements 

BARK number of non-zero elements in the master 
stiffness matrix 

+RHS sum of BARK and the non-zero elements (forces) 
on the right-hand sides of the deflection 
equations 

REDU non-zero elements remaining in +RHS after 
reduction to account for support conditions 

NOTES: 

1. The first page of output is a listing in block tabular 


form of the input values for the X-coordinates. The 
coordinates are listed row-by-row, from left to right, 
in order of the node-point number to which they 
correspond. Index numbers at the top and to the left 
of the block assist the user in identifying the node 
number for each coordinate value. If the node is fully 
restrained by a support in the X-direction, a letter §S 
appears following the value. The second page gives 
Y-coordinate data in a similar format. These data are 
followed by self-explanatory listings of the remaining 
input data. 


The calcuated output begins with block tabular listings 
of deflections at each node point. These deflections 

are relative to the datum coordinate system used to 
define the input coordinates. Tables of calculated 
stresses are then presented for each element, row-by-row. 
The element identification number is shown at the extreme 
left of each row. The stresses are calcalated at the 
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centroid position for each element. 


Next, forces in datum coordinates are listed at each node 
point. These forces are calculated from the computed 
deflections and thus represent a statement of forces 
which must be applied at each node to produce the calcu- 
lated deflections. Generally, a small force will be 
found, even at nodes which were supposed to be unloaded. 
The magnitude of the force is an indication of numerical 
calculation errors. In some cases, oversight errors in 
the input make themselves apparent in these tables. 
Additional checking information is provided by a check 
row which sums the node forces calculated from deflections 
and also sums the moment of these forces about the origin 
of the datum coordinates. 


The last items of output provide information on the 

initial number of non-zero elements contained in the 
simultaneous equations which were solved to find the 
deflections. 


EXAMPLE 1 


This example shows the preparation of input data and provides 
a short problem suitable for testing the operation of the 
program on the user's equipment. Table I shows a listing of 
the input data for the analysis of the structure and loadings 
shown in Figure 4. 


The running time for this example is 1 minute on the IBM 
360 model 40 computer. 
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SHORT ae MUT CASE FUR “ADA 
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Figure 4. 


TABLE I. LISTING OF INPUT CARDS FOR EXAMPLE 1. 


Structure for Example l. 
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Case 2 


EXAMPLE 2 


This example tests the ability of MA2B to find the stresses 
which result from an interference fit between two concentric 
circular cylinders. One-quarter of the assembly was analyzed 
using the finite-element pattern shown in Figure 5. This 
pattern is quite similar to those used for the lug analyses 
reported in Volume I. 


Figure 5 also shows a comparison of the calculated stresses 
with the exact solution. The stresses from MA2B were in 
error by less than 2 percent for this case. The complete 
computer output for this example is shown in Table II. These 
data show a very symmetrical solution. The data also show 
that the interference (i/D = .001) was introduced by speci- 
fying a temperature rise of 134.19F to the inner cylinder. 

No mechanical loads were applied. 


The running time for this example is 10 minutes on the IBM 
-360 model 40 computer. 
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EXAMPLE 3 


This example shows an analysis which is typical of those used 
to find the K,, factors reported in Volume I. The lug consid- 
ered herein is loaded by a force of 1 kip in a direction 45° 
to the axis of symmetry of the lug. In the analysis it is 
assumed that only radial compressive forces exist between the 
lug and the liner. This corresponds to a perfectly greased 
liner. 


Table III presents the entire output fror the computer. The 
output is made more easily understandable by Figure 6, which 
shows the finite-element pattern and the numbering system. 
The node numbezs are shown small; the element numbers, large 
and bold. For clarity, the radial bar elements which connect 
the pin to the liner, and the liner to the lug, are shown 
slightly curved; and the radial boundary positions of the 
liner have been shown with gaps. Two elements, numbers 266 
and 267, are not shown. Their purpose is to provide a path 
for resistance to a tangential force, and thus avoid a singu- 
lar condition in the master stiffness matrix. 


The data in Table III show that a gap was permitted to occur 
between the liner and the lug by specifying a very low modulus 
(only 1 ksi) for all radial bar elements that were in tension. 
The proper assignment of moduli was determined by trial and 
error. 


The tangential stress at the bore of the lug was found by 
extrapolation of plots of tangential stress along radial lines. 
Figure 7 shows these plots. Then the tangential stresses at 
the bore were plotted versus angular position, as shown in 
Figure 8. The magnitude and location of the maximum stress 
can be read from Figure 8. 
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SOURCE FROGRAM 


PROCESSING INFORMATION 


Program MA2B is written in FORTRAN IV (E level) for the IBM 
system 360. It has been run on a model 40 using the Disc 
Operating System Version 3, Level (release 17). Four tape 
drives are required for temporary storage. It can be run with- 
out change on any IBM 360 having 128,000 bytes of storage. It 
may also be run on any computer having a FORTRAN IV compiler 
with only minor changes, providing sufficient storage is 
available. The program consists of a main program plus four 
subroutines. 


An approximate running time for an IBM 360 model 40 computer 
can be found from the following expression: 


Time = (number of nodes) 2/1000 minutes 


FLOW DIAGRAM FOR MA2B 


[Coordinates, Mech. Loads 


—_—___—» Write 13 Mech. Loads 


 eaeaenaons 
5 g | Element Topology 
bal—1. Properties & Temperatures 
+» fq! 
oe | 
| Element Stiffness 
ee, | Master Stiffness | 

' 


Pe pte es aot —-Write 14 Element Data 
——— Write 13 (Mech. & Ther- 

mal Loads) 
Write 10 Master Stiffness 


Add Supports and Forces 
to Master Stiffness 


Deflections 


Stresses 


ee a Temporary Storage 
Read 14 
— Read 10 
Read 13 


Checks | 
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SOURCE PROGRAM DESCRIPTION 


FORTRAN DESCRIPTION 

A matrix cf compatible strain distribution 
for unit element displacements 

AL direction cosines for PQ direction 

AL2 direction cosines for TR direction 

BARK upper part of master stiffness matrix; also 
DBARK and DBAR 

BLK blank symbol on one of input cards 

Cc matrix occurring in Hooke's Law 

CGE an integration over the area of an element, 
used in development of stiffness matrix 

DDSK element stiffness matrix, datum coordinates 

DSK element stiffness matrix, local coordinates 

DQRU element thermal force matrix, datum coordinates 

ECH matrix used for equilibrium checks 

F matrix which relates assumed arbitrary coef- 
ficients in displacement function to the 
displacements 

IBARK identification number for non-zero element 
in master stiffness matrix 

JLAM dimension of ALAM array 

JLAM2 JLAM/2 

MPQRS an array which contains the scheme for build- 
ing the master stiffness matrix 

NBC identification of rows in master stiffness 
matrix corresponding to supports 

NCR@SS humber of supports 

NEL N2 + 1 
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FORTRAN Description 


NUM number of elements in upper half of master 
stiffness matrix 


N2 number of nodes X 2 

QBAR array into which the mechanical loads are read 
at the beginning of the program. The thermal 
loads are subtracted as they are calculated. 
When all elements are processed, QBAR = sum 
of mechanical and thermal loads. 

Q@RU element thermal force matrix, local coordinates 


SCALEF scale factor for master stiffness matrix, 
usually 1.0-5 


STRESS stresses 
Notes: 


1. Several source program symbols were defined in the list 
of input symbols. 


2. All other source symbols are temporary storage. 
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